**************************************************
*        COCHRAN-ARMITAGE TEST FOR TREND         *
**************************************************
* Ref: Alan Agresti (2002), Categorical Data     *
* Analysis (2nd Ed), pp 181-2, John Wiley & Sons *
**************************************************
* (c) Marta Garcia-Granero (september 2005)      *
* biostatistics@terra.es                         *
* downloaded from: http://www.spsstools.net      *
* Feel free to use or modify this code, but      *
* acknowledge the auther and website             *
**************************************************

* Dataset #1: scoring is OK (example provided with SAS at:) 
  (http://www.id.unizh.ch/software/unix/statmath/sas/sasdoc/stat/chap28/sect40.htm)
* (SAS results):
* Statistic (Z)  4.7918
* Asymptotic Test
*  One-sided Pr < Z <.0001
*  Two-sided Pr > |Z| <.0001
* Sample Size = 161

DATA LIST LIST/dose adverse count (3 F8.0).
BEGIN DATA
0 1 26
0 2  6
1 1 26
1 2  7
2 1 23
2 2  9
3 1 18
3 2 14
4 1  9
4 2 23
END DATA.
VALUE LABEL adverse 1'No' 2'Yes'.
WEIGHT  BY count .

* Using SPSS linear by linear test (close enought) *.
CROSSTABS
  /TABLES=dose BY adverse
  /FORMAT= AVALUE TABLES
  /CELLS= COUNT ROW
  /STATISTIC=CHISQ.

* Cochran-Armitage test with MATRIX *.

* This code assumes the dependent variable is coded 0/1 (No/Yes) *.
TEMPORARY.
COMPUTE adverse=adverse-1.
AGGREGATE
  /OUTFILE = 'c:\temp\aggdata.sav'
  /BREAK = dose
  /pi = MEAN(adverse)
  /ni = N.

MATRIX.
PRINT/TITLE='Cochran-Armitage Test for Trend (Agresti 2002)'.
GET pi /VAR=pi /FILE='c:\temp\aggdata.sav'.
GET ni /VAR=ni /FILE='c:\temp\aggdata.sav'.
* Scoring (change if required) *.
COMPUTE xi={0; 1; 2; 3; 4}.
* Reports *.
PRINT {xi,ni&*pi,ni}
 /FORMAT='F8.1'
 /TITLE='Statistics'
 /CLABELS='Score i','Count i','N i'.
PRINT MSUM(ni)
 /FORMAT='F8.0'
 /TITLE='Total sample size'.
* Calculations *.
COMPUTE p=MSUM(ni&*pi)/MSUM(ni).
COMPUTE x=MSUM(xi&*ni)/MSUM(ni).
COMPUTE num=(xi-x)&*ni&*pi.
COMPUTE den=p*(1-p)*ni&*((xi-x)&**2).
COMPUTE z=MSUM(num)/SQRT(MSUM(den)).
COMPUTE zsig=2*(1-CDFNORM(ABS(z))).
* Final report *.
PRINT {z,z**2,zsig,zsig/2}
 /FORMAT='F8.3'
 /CLABELS='Z','Z˛','2-sig','1-sig'
 /TITLE='Statistics & Significance'.
COMPUTE observed={(ni-ni&*pi),ni&*pi}.
COMPUTE expected=RSUM(observed)*CSUM(observed)/MSUM(ni).
COMPUTE chisq=MSUM((observed-expected)&**2/expected).
COMPUTE chisqp=1-CHICDF(chisq,NROW(observed)-1).
COMPUTE reschi=chisq-z**2.
COMPUTE reschip=1-CHICDF(reschi,NROW(observed)-2).
PRINT {chisq,chisqp;reschi,reschip}
 /FORMAT='F8.3'
 /RLABELS='Pearson','Deviat.'
 /CLABELS='Chi˛','Sig.'
 /TITLE='Pearson-Chi˛ & Deviation from trend (non-linearity)'.
END MATRIX.

* As you can't specify scores for CROSSTABS, the key is recoding
* the predictor to asign the scores you want. Suppose your data
* don't reflect the scoring you want to give *.

* Dataset #2 (Agresti 2002) *.
DATA LIST LIST /alcohol malform count (3 F8.0).
BEGIN DATA
0 1    48
0 0 17066
1 1    38
1 0 14464
2 1     5
2 0   788
3 1     1
3 0   126
4 1     1
4 0    37
END DATA.

VALUE LABEL alcohol
 0 ' 0'
 1 '<1'
 2 '1-2'
 3 '3-5'
 4 '>6'.
VALUE LABEL malform 0'Absent' 1'Present'.
WEIGHT  BY count .

* Agresti gives the following scores: 0, 0.5, 1.5, 4, 7 *.

* The following will not give the results expected (Z˛=6.57 p=0.010): *.

CROSSTABS
  /TABLES=alcohol BY malform
  /FORMAT= AVALUE TABLES
  /CELLS= COUNT ROW
  /STATISTIC=CHISQ. /* LxL test is 1.828 (p=0.176)*. 

* Recoding and re-running *.
TEMPORARY.
RECODE alcohol (0=0) (1=0.5) (2=1.5) (3=4) (4=7) INTO alcscore .
VARIABLE LABELS alcscore 'Scores for alcohol'.
* Using SPSS linear by linear test (close enough) *.
CROSSTABS
  /TABLES=alcscore BY malform
  /FORMAT= AVALUE TABLES
  /CELLS= COUNT ROW
  /STATISTIC=CHISQ. /* LxL is now 6.570 *.

* Cochran-Armitage test with MATRIX *.

* This time, MALFORM takes the correct values (0/1) *.
AGGREGATE
  /OUTFILE='c:\temp\aggdata.sav'
  /BREAK=alcohol
  /pi = MEAN(malform)
  /ni=N.

MATRIX.
PRINT/TITLE='Cochran-Armitage Test for Trend (Agresti 2002)'.
GET pi /VAR=pi /FILE='c:\temp\aggdata.sav'.
GET ni /VAR=ni /FILE='c:\temp\aggdata.sav'.
* Scoring (change if required) *.
COMPUTE xi={0; 0.5; 1.5; 4; 7}.
* Reports *.
PRINT {xi,ni&*pi,ni}
 /FORMAT='F8.1'
 /TITLE='Statistics'
 /CLABELS='Score i','Count i','N i'.
PRINT MSUM(ni)
 /FORMAT='F8.0'
 /TITLE='Total sample size'.
* Calculations *.
COMPUTE p=MSUM(ni&*pi)/MSUM(ni).
COMPUTE x=MSUM(xi&*ni)/MSUM(ni).
COMPUTE num=(xi-x)&*ni&*pi.
COMPUTE den=p*(1-p)*ni&*((xi-x)&**2).
COMPUTE z=MSUM(num)/SQRT(MSUM(den)).
COMPUTE zsig=2*(1-CDFNORM(ABS(z))).
* Final report *.
PRINT {z,z**2,zsig,zsig/2}
 /FORMAT='F8.3'
 /CLABELS='Z','Z˛','2-sig','1-sig'
 /TITLE='Statistics & Significance'.
COMPUTE observed={(ni-ni&*pi),ni&*pi}.
COMPUTE expected=RSUM(observed)*CSUM(observed)/MSUM(ni).
COMPUTE chisq=MSUM((observed-expected)&**2/expected).
COMPUTE chisqp=1-CHICDF(chisq,NROW(observed)-1).
COMPUTE reschi=chisq-z**2.
COMPUTE reschip=1-CHICDF(reschi,NROW(observed)-2).
PRINT {chisq,chisqp;reschi,reschip}
 /FORMAT='F8.3'
 /RLABELS='Pearson','Deviat.'
 /CLABELS='Chi˛','Sig.'
 /TITLE='Pearson-Chi˛ & Deviation from trend (non-linearity)'.
END MATRIX.